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Abstract. Hydrodynamical simulations are presented of a pulsar wind emitted by a supersonically moving pulsar. 
The pulsar moves through the interstellar medium or, in the more interesting case, through the supernova remnant 
created at its birth event. In both cases there exists a three-fold structure consisting of the wind termination shock, 
contact discontinuity and a bow shock bounding the pulsar wind nebula. Using hydrodynamical simulations we 
study the behaviour of the pulsar wind nebula inside a supernova remnant, and in particular the interaction 
with the outer shell of swept up interstellar matter and the blast wave surrounding the remnant. This interaction 
occurs when the pulsar breaks out of the supernova remnant. We assume the remnant is in the Sedov stage of its 
evolution. Just before break-through, the Mach number associated with the pulsar motion equals M pBI = 7/V5, 
independent of the supernova explosion energy and pulsar velocity. The bow shock structure is shown to survive 
this break-through event. 

Key words. Pulsars - Supernova Remnants - Shocks - Hydrodynamics 



1. Introduction 

A supernova explosion of a massive star will result in an 
expanding supernova remnant (SNR). In some cases the 
fossil of the progenitor star is a pulsar moving at high 
velocity. Even though the precise physical mechanism re- 
sponsible for imparting a large kick velocity to single ra- 
dio pulsars at birth has not been identified, observations 
of the pulsar distribution with respect to the mid-plane of 
the galaxy indicate that they are born with a velocity in 
the range V psr ~ 100 - 1000 km/s (Harrison et al. 1992; 
Lyne & Lorimer 1994). A similar range of values is ob- 
tained from a sample of SNR-pulsar associations (Frail et 
al. 1994). 

The expansion of a SNR is decelerated due to mass- 
loading by swept up interstellar medium (ISM) or by ma- 
terial from a progenitor wind. As the pulsar moves with 
a constant velocity it will ultimately break through the 
SNR shell. Two observed systems are often presented as 
an illustration of this scenario: 



CTB80: in this supernova remnant the pulsar PSR 
B1951+32 is located (in projection) just inside the outer 
edge of the remnant. The spectral index of the synchrotron 
emission in the vicinity of the pulsar system indicates that 



there is a plerionic nebula around the pulsar, see for ex- 
ample Strom (1987) and Strom & Stappers (2000). 

G5. 4-1.2: in this case the pulsar is located well outside 
the supernova remnant. At radio frequencies an emission 
bridge appears to connect the pulsar B 1757-24 an d the 
associated pulsar wind nebula (PWN) with the supernova 
remnant (Frail & Kulkarni 1991), suggesting a physical 
association between the supernova remnant and the pul- 
sar. It should be pointed out, however, that a new upper 
limit on the proper motion of Bl 757-24 (Gaensler & Frail 
2000) puts the component of the pulsar velocity in the 
plane of the sky at Vj_ psr < 600 km/s for an assumed 
distance of 5 kpc. This leads to a discrepancy between 
the characteristic pulsar age, obtained from its spin pe- 
riod derivative (P/2P ~ 16 kyr), and the dynamical age 
obtained from the offset distance i? psr from the center of 
G5.4-1-2 (R psr /V psr > 39 kyr). 

Both systems are clearly brightened at radio wave- 
lengths near the position of the pulsar, and it has been 
suggested that the associated pulsar wind is rejuvenating 
the radio emission from the SNR shell by the injection of 
fresh relativistic electrons (Shull et al. 1989). In this pa- 
per we will investigate the hydrodynamical aspects of the 
interaction between a pulsar wind and a SNR shell. 
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Most pulsars have a lifetime (10 6 — 10 7 yr) which is 
much larger than the age < 10 4 yr of a SNR in the Sedov 
phase. Therefore pulsars will remain visible long after the 
associated SNR has dissolved into the interstellar medium. 
The pulsar then moves as an isolated pulsar through the 
interstellar medium, and can form a pulsar wind nebula 
bow shock system. A typical example of such a system is 
the Guitar Nebula around PSR B2224+65 which has been 
detected both in Ha (Cordes et al. 1993) and in X-rays 
(Romani et al. 1997), but which has no associated SNR. 

In this paper we consider the case where the pulsar's 
kick velocity is sufficiently high so that it leaves the super- 
nova remnant while it is still in the Sedov stage. We de- 
scribe three different stages in the evolution of the pulsar- 
SNR system: (1) the stage where the PWN/bow shock 
resides inside the SNR, (2) the PWN/bow shock breaking 
through the shell of the SNR and (3) the stage where the 
PWN/bow shock moves through the ISM. 



i PDS = 2.9 x 10 4 Et{ 17 n- 9/17 yr. (1) 

Here -E51 is the explosion energy E snl - in units of 10 51 ergs 
and hq denotes the hydrogen number density in the ISM. 

We will describe the physics of a pulsar wind interac- 
tion with the shell of a SNR in the Sedov stage. The results 
presented below only apply for certain range of values for 
the pulsar velocity V^ sr . Equating the distance traveled by 
the pulsar, 

Rpsr — ^psr^ > 

with the radius for a SNR embedded in a homogeneous 
interstellar medium of density pi sm in the Sedov stage, 

/ f \ 1 / 5 

i? snr ^1.15p^ < 2 / 5 , (2) 



2. Physics of a PWN bow shock inside a SNR 

2.1. Dynamics of the pulsar/SNR system 

In rapidly rotating (young or recycled) pulsars, it is be- 
lieved that a pulsar wind carries away most of the spin- 
down luminosity, 



one gets the crossing time for the pulsar: 



t cr = 1.27 



which is 



Pisml^sr 



1/3 



(3) 



l = mn , 

of a pulsar with rotation period P — 2ir/Q and moment of 
inertia I. This relativistic wind is presumably generated 
in the pulsar magnetosphere, and accelerates electrons, 
positrons and possibly nuclei to ultra-relativistic speeds. 

The pulsar wind blows a bubble or pulsar wind nebula 
(PWN) into the surrounding medium. This PWN is ini- 
tially located well within the interior of the SNR created 
at the birth of the neutron star. During the free expansion 
stage of the SNR evolution the typical expansion speed of 
the stellar ejecta as determined by the mechanical energy 
E sm released in the explosion and the ejecta mass M C j, 



'10 E s . 



3 M, 



10 000 km/s , 



1.4 x 10 4 EU 3 V 10 U 3 n 



1/3 

y- 



(4) 



Here V1000 is the velocity of the pulsar in units of 1,000 
km/sec, and n the hydrogen number density of the ISM 
in units of cm -3 . The requirement t cr < tpos yields the 
minimum velocity a pulsar needs in order to break out of 
the SNR while the latter is still in the Sedov stage: 



V psr > 650 n 



2/17^,1/17 



51 



km/s. 



(5) 



Although this is a rather high value, it is still in the range 
of values observed by Harrison et al. (1992). 

One can use the Rankine-Hugoniot relations to deter- 
mine the pressure just behind the Sedov- Taylor blast wave 
bounding the SNR (assuming a strong shock and a gas 
with specific heat ratio 7 = 5/3) 



is generally much larger than the kick velocity of the pul- 
sar. As a result the PWN is located relatively close to 
the center of the SNR at this stage. In this regime, the re- 
sults of our earlier, spherically symmetric simulations (van 
der Swaluw et al. 2001) are expected to remain approxi- 
mately valid. Only when the SNR expansion slows down 
as it enters the Sedov stage after some ~ 500 — 1 000 yr, 
a situation is possible where the pulsar position becomes 
strongly excentric with respect to the SNR. 

The Sedov stage of SNR expansion lasts until internal 
(radiative) cooling becomes important. The SNR then en- 
ters the so-called pressure-driven snowplow (PDS) stage. 
The relevant transition time is calculated by Blondin et 
al. (1998): 



-Psh = 



3- o- V 2 

a /Asm v am 



(6) 



where 



V m r = 



dR 

snr 

dt 



2 ± <-snr 

5 — r~ 



is the SNR expansion speed. Using this expression to- 
gether with expression (|3|) for the crossing time and the 
Sedov solution (^|) , one can derive the following equations 
valid at the moment of break-through. The speed of the 
pulsar is related to the SNR expansion speed by 



Vry. 



5 V 

o * SI 



(7) 
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Fig. 1. Configuration of the pulsar wind nebula moving 
through a uniform medium as seen in the rest frame of 
the pulsar. 



while the material in the shell behind the SNR blast wave 
moves with a velocity 



3 V 

A * SI 



(8) 



This corresponds to a relative speed between pulsar and 
post-shock material equal to 



V Ie \ = V psr — Vgh. — T V* 



4 'snr • 



(9) 



The density in the shell is roughly the post-shock density 
for a strong shock, p s h = 4p ism , so the Mach number M pS r 
of the pulsar motion through the shell material satisfies 



Mr 



v / 7 P sh /4p if 



7_ 

71 



3.13 



(10) 



Since the proper motion of the pulsar is supersonic with 
respect to the surrounding medium the outer rim of the 
PWN will become a-spherical and a bow shock will form. 

2.2. Pulsar Wind 

A pulsar wind is believed to consist of an ultra-relativistic, 
cold flow with a large bulk Lorentz factor, e.g. T w > 10 3 in 
the case of the Crab (Gallant et al. 2002). The cold wind 
is terminated by a termination shock which thermalizes 
the flow, leading to a relativistically hot downstream state 
with sound speed s ~ c/V3. Following Kennel & Coroniti 
(1984), the luminosity of the pulsar wind at the wind ter- 
mination shock, Rt s is given by a combination of particle 
and magnetic luminosity: 



L = 47rr^n w i^ 8 mc 3 (l + a) 



(11) 



where n w is the proper density in the wind, m the mean 
mass per particle and a is the ratio of Poynting flux to 
plasma kinetic energy flux. In this paper we consider an 
unmagnetised pulsar wind, i.e. a = 0. The typical pres- 
sure behind the ultra-relativistic termination shock is (e.g. 
Kennel & Coroniti 1984) 



2 p2 2 



L 



6itR£ 



(12) 



The first equality is approximate because of deviations of 
sphericity of the pulsar wind region, induced by the proper 
motion of the pulsar. The shocked pulsar wind material is 
separated from material that has gone through the bow 
shock by a contact discontinuity. Because of the high inter- 
nal sound speed, both in the pulsar wind material and in 
the material that has passed through the bow shock, and 
because of the small size of the region between the termi- 
nation shock and bow shock, this region can be considered 
to be isobaric to lowest approximation. Using the Rankinc- 
Hugoniot jump conditions at the stagnation point at the 
head of the bow shock surrounding the PWN, the pressure 
equals: 



P bs = (57W^ sr - l)P sh /4 = 9 p lsm ^ r . 



(13) 



After the pulsar has broken through the shell the pulsar 
wind is completely confined by the ram pressure of the cold 
ISM and the stagnation-point pressure, using the Rankine- 
Hugoniot relations again, drops to 



Pbs — 4Pism^C 



psr 



16 Pism V snr j 



(14) 



a pressure reduction by roughly 50% as the pulsar leaves 
the SNR. The fact that the region between termination 
shock and bow shock is almost isobaric implies 



ft: 



This condition, together with Eqn. (|lj), determines the 
stand-off distance of the pulsar wind termination shock as 



Rts = V 



L 



6irp ism V 2 sl .c 



1/2 



(15) 



where the numerical factor rj is determined by Eqns ( [13] ) or 
(pT|). This parameter takes the value rj = 5/6 ~ 0.83 when 
the pulsar is still just inside the SNR, and 77 = 2/v3 ~ 
1.15 when the pulsar moves through the ISM. The radius 
Rts is also the typical stand-off distance of the bow shock, 
which is always close to the termination shock at the head 
of the pulsar wind nebula. 

These expressions allow us to calculate the relative size 
of the pulsar wind to the supernova remnant at the mo- 
ment of break-through. From the expression (|j|) for the 
crossing time one has 
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Table 1: Simulation parameters 



Parameter 


Simulation 


Pulsar wind luminosity Lq (erg/s) 


1.3 x 10 34 


Mass density po (g/cm 3 ) 


1.0 x 10" 24 


Terminal Velocity (cm/sec) 


3.0 x 10 9 


Velocity Medium V re \ (km/sec) 


500.0 


Mach Number Pulsar Ai psT 


3.13 


AT "U £ 'J 11/' 1 * j ■ \ 

JN umber or grid cells (m r-directionj 


300 


Number of grid cells (in z-direction) 


180 


Grid size (pc) (in r-direction) 


0.1 


Grid size (pc) (in z-direction) 


0.06 



Rsnri^cr) — ^psr^cr — 13.6 ^lOOO n P*" ' 

The termination shock radius is of order 

R ts ~ 2.8 x 10~ 3 r) L% 2 n a 1/2 V w l pc , 

here L 36 = L/(10 36 erg/s). Note that the size of the SNR 
shell is much larger then the size of the PWN. For this rea- 
son we will neglect the curvature of the SNR blast wave 
and perform a hydrodynamical simulation where the pul- 
sar moves with a Mach number of A^psr = 7/V5 ~ 3.13 
through the post-shock flow of a strong plane-parallel 
shock, ultimately crossing this shock into the unshocked 
medium. 




0.00 

0.00 0.02 0.04 0.06 0.08 0.10 
Radius (in parsec) 

Fig. 2. Comparison between the numerical result for the 
bow shock with a low Mach Number (M = « 3.13) 

with the profile given by Wilkin (1996) for M. 3> 1 (solid 
line). The grey scale represents the pressure. 

3.2. Initialising a pulsar wind 

The current version of the hydrodynamics code uses one 
fluid with a single equation of state in the computa- 
tional domain. For simplicity, and accuracy in the non- 
relativistic part of the flow, we use an ideal fluid with spe- 
cific heat ratio 7 = 5/3. The pulsar wind is modeled by 
continuously depositing thermal energy at a constant rate 
L (the spindown luminosity) in a small volume, together 
with an associated mass injection M pw . The hydrodynam- 
ics code itself then develops a thermally driven wind with 
terminal velocity which is subsequently randomised at 
a termination shock. In order to increase the computation- 
ally efficiency, we take the mechanical luminosity L and 
mass deposition rate M pw such that the terminal velocity 
of the wind as determined from these two parameters is 
much larger than the other velocities of interest: 



3. Numerical simulations of the PWN bow shock 

3.1. Simulation Method 

We simulate a pulsar wind using the Versatile Advection 
Code [] (VAC), a general -purpose software package devel- 
oped initially by G. Toth at the Astronomical Institute 
in Utrecht (Toth 1996; Toth & Odstrcil 1996). The con- 
figuration is depicted in figure 1, showing both shocks 
which are of interest: the pulsar wind termination shock 
and the bow shock bounding the PWN. This system is 
assumed to be axially symmetric around the direction 
of motion of the pulsar. Out of the several choices for 
discretizing the equations of hydrodynamics in conserva- 
tive form available in VAC, we use a shock-capturing, 
Total- Variation-Diminishing Lax-Friedrich scheme (Toth 
& Odstrcil, 1996). Our simulations have been performed 
in two dimensions, using a cylindrical coordinate system. 



= J2L/M ] 



pw ■ 



This choice will result in the correct global behaviour of 
the PWN. This method is similar to the method as de- 
scribed in van der Swaluw et al. (2001). We employ a 
non-uniform grid with highest resolution near the pulsar 
in order to fully resolve the pulsar wind region. The non- 
uniform grid is chosen in such a way that the difference in 
size between adjacent grid cells is always less then 5 % . 
The radius of the termination shock, given by Eqn. (15) 



in the relativistic case, is replaced by its non-relativistic 
equivalent, 



Rt, 



L 



1/2 



1 See 



http://www.phys.uu.nl/~toth/ 



and will have roughly the correct value when ~ c. 
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Fig. 3. Gray-scale representation of the density distribu- Fig. 5. Gray-scale representation of the pressure distribu- 
tion of a PWN and bow shock with the parameters as tion of a PWN bow shock with the same parameters as in 
denoted in table 1. Figure 3. 




-0.030 -0.020 -0.010 0.000 0.010 0.020 0.030 
Z (in parsec) 

Fig. 4. Density profile of the PWN and bow shock along 
the z-axis of Figure 3. The location of the Pulsar, forward 
and backward termination shocks and of the bow shock are 
indicated. The pulsar wind region is shaded light gray, the 
region containing shocked interstellar material is shaded 
a darker gray. 
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Fig. 6. Pressure profile of a PWN and bow shock along 
the z-axis of figure 5. Shading is as in Figure 4. 

ables at time t n , denoted by at grid point i, with 
their values at the previous time t n —i. We then calculate 
the residual Res defined as 



3.3. Steady PWN bow shock in a uniform medium 

Our calculations are performed in the pulsar rest frame. 
The pulsar wind nebula is allowed to evolve in a uniform 
medium, moving at a constant speed V psr at large dis- 
tances from the pulsar. This medium represents the in- 
terior of the supernova remnant (shocked ISM) close to 
the blast wave. The velocity V psr is supersonic with re- 
spect to the internal sound speed of the medium so that a 
bow shock develops around the PWN. We let the hydro- 
dynamics code evolve the system until the large-scale flow 
is steady. 

This simulation has been performed with parameters 
as denoted in table 1. In order to determine when the 
system is steady, we employ the prescription of Toth et 
al.(1998). This prescription compares all iV var flow vari- 



Res 



\ 



iV v; 



JV var 

E 

a=l 



n+1 



U n -\ 



(16) 



and halt the calculation when Res has a value less than a 
predetermined critical value, typically Res = 1CP 4 . 

Wilkin (1996) has given an analytical equation for the 
geometry of a stellar wind bow shock. His solution, for the 
distance r to the wind source in terms of the polar angle 
9 with respect to the symmetry axis, reads: 



(17) 



Here Rq « i?bs is the stand-off distance of the bow 
shock on the symmetry axis (9 = 0). We compare our 
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Fig. 7. Pressure distribution of a PWN bow shock with 
parameters as denoted in table 1. Here the PWN is inter- 
acting with the shell of the SNR, modeled as a plane- 
parallel shock. The profile given by Wilkin (1996) for 
M. S> 1 (solid line) agrees well with the bow shock profile 
in the ISM. 

morphology with Wilkin's result, where we equate Ro to 
the stand-off distance of the bow shock in the simula- 
tions. This is depicted in hgure 2. The asymptotic cone 
of Wilkin's solution is seen to be significantly narrower 
than that of the bow shock in the simulation. This is be- 
cause his solution balances only the ram pressures of the 
wind and the ambient medium, i.e. it corresponds to the 
limit M psr » 1, while in our case the Mach number is 
moderate: M psr « 3.13. 

The figures 3 and 4 show density profiles of the PWN 
bow shock of a pulsar moving through a uniform medium. 
One can see the contact discontinuity, midway between 
the termination shock and the bow shock, separating 
the shocked pulsar wind material and the much denser 
shocked ISM. The synchrotron emission of the plcrionic 
PWN is expected to come from the shocked pulsar wind 
material, whereas the material swept up by the bow shock 
can show up as Ha emission. 

Figure 5 shows the pressure distribution and figure 6 
shows the pressure profile along the symmetry axis. One 
can see the pulsar wind, originating at the pulsar position 
z = 0, and its termination shock atzK —0.005 pc ahead 
of the pulsar in the direction of motion, and at z ~ 0.017 
pc behind the pulsar. The bow shock bounding the PWN 
is located at z ~ —0.012 pc. The region between the pul- 
sar wind termination shock and the bow shock is almost 
isobaric. As shown by van der Swaluw et al. (2001), this 
is also the case for a PWN around a stationary pulsar 
located at the center of the SNR. 

3.4. Interaction of the PWN with a shock 

In this section we present results of the break-through of 
the PWN bow shock when it crosses the shell of its super- 
nova remnant. The simulation is performed once more in 




0.00 I , , , i , , , i , , , i , , , i , , , I 

0.00 0.02 0.04 0.06 0.08 0.10 
R (in parsec) 

Fig. 8. Density distribution of the PWN bow shock-SNR 
interaction corrsponding to the pressure distribution of 
Figure 7 

the rest frame of the pulsar. We initialise a steady-state 
configuration of the PWN bow shock as described above, 
with the same parameters as denoted in table 1. Next we 
use the Rankine-Hugoniot relations to initiate a strong 
shock front (with Mach number M. ~ 100) near the bot- 
tom (low z) end of the grid, with which the pulsar wind 
bow shock then catches up in accordance with equations 
(7)-(10). As stated in section 2, the PWN bow shock is 
much smaller than the radius of the SNR, so we can safely 
approximate the SNR blast wave as a plane-parallel strong 
shock. 

At the end of the simulation, when the strong shock 
is almost at the upper boundary of the grid, numerical 
instabilities arise. Therefore we stop the simulation after 
the configuration as shown in the figures 7-9, when the 
influence of the numerical instabilities are not influencing 
the solution too strongly. 

Figures 7 and 8 show the system, some time after the 
SNR shock has passed the head of the bow shock. Figure 
7 shows the pressure distribution and figure 8 the density 
distribution. Figure 9, which compares the density pro- 
file along the z-axis before and after the passage of the 
supernova blast wave, shows that the pulsar wind neb- 
ula expands, roughly by a factor 1.5, after it leaves the 
SNR. This reflects the reduction in the confining (ram- 
)pressure, as calculated in Section 2. Figure 7 also shows 
that the bow shock shape in the ISM now closely fits 
the analytical result of Wilkin (1996). Once the pulsar 
wind bow shock has crossed the supernova blast wave, its 
Mach number with respect to the interstellar medium is 
much larger than unity. This agreement is all the more 
remarkable given that strictly speaking, the assumptions 
under which Wilkin's solution is derived (geometrical thin- 
ness and complete mixing of the post-shock fluids) are 
not fulfilled in pulsar bow-shock nebulae (Bucciantini & 
Bandiera 2001). 

During the interaction between the pulsar wind and 
the shell of the SNR, the PWN bow shock and the SNR 
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Fig. 9. Density profile of the PWN bow shock along the 
z-axis. In this figure the pulsar's position corresponds to 
Z = 0. The solid curve denotes the density profile before 
interacting with the SNR shock. The dashed curve denotes 
the density profile well after the SNR shock has passed the 
head of the bow shock. The dashed line and the dash-dot 
line give the position of the bow shock around the PWN 
before and after the passage of the supernova blast wave. 
The PWN expands roughly by a factor 1.5. 



blast wave intersect. This intersection produces an addi- 
tional pressure gradient which results in an accumulation 
of mass. The resulting pressure and density enhancements 
can be seen as dark spots at the region of intersection in 
the figures 7 and 8. When the bow shock moves through 
the shell of the remnant it encounters the unshocked ISM. 
The ambient density is reduced by a factor 4, which results 
in a similar density reduction behind the bow shock. 



4. Conclusions 

We have considered the case of a pulsar wind breaking 
through the shell of a SNR in the Sedov- Taylor stage. Only 
high- velocity pulsars reach the edge of the SNR while the 
latter is still in the Sedov stage of its evolution. At the 
moment of break-through, the ratio of the pulsar velocity 
and SNR expansion speed is V psl /V snr = 5/2, and the 
Mach number associated with the pulsar motion equals 
■M. pST = 7/v5- These conclusions are independent of the 
explosion energy E snr or the pulsar speed V psr . 

Our simulations show that the break-through of the 
PWN does not lead to a significant disruption. The re- 
duction of stagnation pressure by about 50% leads to a 
moderate expansion of the PWN where its radius increases 
by a factor ^1.5. The latter result was also obtained an- 
alytically in Eqns. (13)-(15). 

There is good agreement between our numerical results 
and analytical estimates, based on pressure balance argu- 
ments, for the size of the bow shock surrounding the PWN. 
The only clear indication of the interaction between the 
PWN bow shock and the SNR (Sedov- Taylor) blast wave 



is a density and pressure enhancement at the intersection 
of these two shocks. 

In a subsequent paper we will consider the effects of 
the energetic particles which were injected by the pulsar 
wind into the surroundings; some preliminary results of 
this investigation can be found in van der Swaluw et al. 
(2002). 
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